Week 3: Reaction-diffusion systems

Week 3: Reaction-diffusion systems#

How the cheetah got his spots#

Why aren’t they striped?

Fig: from ChatGPT

Patterns#

And it’s not just a cheetah’s spots… there are patterns all around us.

We see similarities in systems that seem to share nothing in common.

Pattern formation is explored in many fields: developmental biology, cell biology and synthetic biology, to physics, mathematics and computer science.

Complex Systems Thinkers seek to use mathematical methods to:

model and understand the underlying mechanisms that generate the patterns; and
classify and quantify the resulting pattern.

The goal is not to focus on system-specific details, but to uncover general principles and shared dynamics that transcend disciplinary boundaries.

What exactly is a pattern?#

We sort of know intuitively. Humans have an uncanny ability to smooth detail in order to see patterns.

Neither bland uniformity nor disorderly chaos but something in between with some notion of ‘repetition’.

A simple type of pattern is a repeating structure in space, like wallpaper tessellations or symmetries (shifting one’s view by one repeat length results in seeing the same thing). There are also repeating patterns in time, like the seasons.

Alternatively, you can consider something like a dress pattern (prototype), which allows for the dress to be copied.

Possible definitions:

A regularity or structure observed in a system’s behaviour or configuration that arises from the interactions among its parts.
An emergent property of a system that is visible at the macro level but not attributable to any individual component alone.
A compressible description of a system—i.e. a configuration that allows the system to be described more succinctly than by listing each part independently.
A set of statistically or structurally consistent relationships across space, time, or scale within a system or across similar systems.
A form of order that reflects symmetry, repetition, correlation, or predictability, often despite underlying randomness or complexity.

How do patterns form?#

How patterns form (pattern formation) and how they are recognised by the brain (pattern recognition) are important areas of study in Complex Systems.

We have already seen fractals, which show self-similar structure across scales (e.g. in plant branching or river networks).

But there are many other processes in the natural world that form patterns.

Rayleigh–Bénard convection: where heat flow generates stable, symmetric roll patterns in fluids.

The pattern that forms in a hot coffee on a cold morning might look like a bit like this garlic but they are very different mechanisms. The semi-regular pattern of garlic cloves is the result of a growth process producing a such as honeycomb-like structures (Voronoi tesselation) in biological tissues or crystals.

Fig from: @alizauf (on X)

Elastic instabilities: are thought to underlie folding patterns in the cerebral cortex of higher animals

Travelling waves, fronts, oscillations, and clustering: as seen in chemical reactions, neural activity, or population dynamics

from IPython.display import YouTubeVideo
display(YouTubeVideo("jRQAndvF4sM", width=400))
# https://www.youtube.com/watch?v=jRQAndvF4sM

Diffusion-based mechanisms: including reaction–diffusion systems that generate spots, stripes, and spirals.

This is what we will focus on and in particular, one of the most celebrated mechanisms to explain self-organising spatial structures… Turing reaction-diffusion instability as a source of symmetry breaking.

These patterns arise not from a single process in isolation, but from the interplay of multiple mechanisms operating across different spatial and temporal scales.

In other words, they emerge through what we call spatio-temporal dynamics — the evolving behaviour of a system over time and across space. This captures both the when and the where of system change. Note that “space” may refer to physical space (e.g. tissue, fluid) or state space (e.g. neural activity, chemical concentration).

Understanding spatio-temporal dynamics is essential in Complex Systems science because many hallmark behaviours, such as pattern formation, synchronisation, and self-organisation—unfold across both space and time. These behaviours cannot be explained by static snapshots or purely local interactions alone.

Spatio-temporal dynamics are typically modelled using partial differential equations (PDEs), which describe continuous changes across fields. This distinction will become important next week when we explore cellular automata, which take a discrete approach to modelling space, time, and state.

Alan Mathison Turing#

Turing is widely recognised as one of the most influential thinkers of the 20th century.

He is known for groundbreaking contributions across multiple domains:

Codebreaking hero: Played a pivotal role in decrypting the German Navy’s Enigma code during the Second World War. His work at Bletchley Park is credited with significantly shortening the war and saving millions of lives.
Founding figure of computer science: Often called the father of computer science, Turing proposed the idea of a universal computing machine—a theoretical blueprint for the digital computer capable of storing and executing programs. His ideas laid the foundation for modern computing.
Visionary in artificial intelligence: Turing also pioneered thinking about machine intelligence. He proposed what became known as the Turing Test—a way of evaluating whether a machine can exhibit human-like intelligence by engaging in a conversation indistinguishable from that of a human. This was famously used in Ridley Scott’s 1982 movie Bladerunner to identify non-human ‘replicants’.
Persecuted for his sexuality: Despite his immense contributions, Turing was subjected to state-sanctioned punishment for homosexual acts, which were criminalised at the time. He accepted chemical castration over imprisonment, and died by suicide in 1954. The British government issued a formal apology in 2009 for its shameful treatment of him.

The chemical basis of morphogenesis#

Turing was ahead of his time in many ways.

Although less well-known for it, he made a profound contribution to a fundamental biological question: how do regular patterns like stripes, spots, and spirals, form in living organisms?

In his 1952 paper “The Chemical Basis of Morphogenesis”, he was the first to propose a mathematical and chemical explanation for this process: that simple chemical interactions—specifically reaction and diffusion—could give rise to complex, self-organising spatial patterns during early development.

This idea was revolutionary. Almost no one was asking these questions at the time; his paper cites just six relevant works! Published in the year prior to his death, this work has been cited over ~~16000 (in 2023)~~ 18000 (in 2025) times and is now recognised as a foundation for the study of biological pattern formation, and a landmark in both mathematical biology and Complex Systems.

Turing, Alan M. The Chemical Basis of Morphogenesis. 1952. Philosophical Transactions of the Royal Society; 237: 37-72

Reaction-diffusion systems#

Turing proposed that complex spatial structures may arise from the interaction of two types of chemical substances, ‘morphogens’ diffusing through tissue.

‘Morphogenesis’, from the Greek morphê meaning shape/form and genesis meaning creation (literally “the generation of form”), is a general mechanism producing a diverse array of patterns.

Morphogens are the ‘shape formers’. They are different chemical signals/substances that react together and diffuse through a plant or animal tissue, one activating and one deactivating growth, to set up patterns of development.

Turing was deliberately vague about what these morphogens were. They could be chemicals/molecules/proteins/hormones/genes etc, operating in tissue or cells… It was the 50s so lots wasn’t known yet (and Turing wasn’t trying to be a biologist).

It is the interaction of two morphogens: an “activator” and an “inhibitor” with competing interests that allows for the autonomous generation of spatial patterns during development via a type of symmetry breaking.

This ‘dappled’ pattern in was created by Turing ‘in a few hours by a manual computation’:

Activator: a slow-diffusing chemical that promotes local change—such as stimulating pigment production in a region of tissue. It acts as a positive feedback signal:
- Encourages the cell to produce pigment (or some morphological change)
- Promotes its own production
- Triggers production of its own inhibitor
Inhibitor: a fast-diffusing chemical that suppresses the effect of the activator. It spreads over a wider area and provides long-range negative feedback:
- Inhibits the production of the activator
- By doing so, it ultimately limits its own production as well

Together, this creates a dynamic where local activation and long-range inhibition interact to produce stable, spatially periodic patterns.

Reaction#

We use \(A\) and \(B\) to denote the activator and inhibitor respectively and will explore the specific example of the following chemical reactions:

\[ A+2B\rightarrow 3B \]

i.e. a single species of \(A\) will interact with two species of \(B\) to form three species of \(B\) (sort of like predator-prey system)

This is part of the Gray-Scott model.

For this interaction to happen (with some probability) the particles first have to come into contact with each other, which happens via diffusion.

Diffusion#

But first… Brownian motion (aka Wiener processes) was named for Robert Brown (1827) who noticed that pollen grains suspended in water moved in a move in a jittery, random fashion.

The underlying cause remained a mystery until Albert Einstein (1905) explained it using the kinetic theory of gases, providing a crucial bridge between the microscopic world of atoms and the macroscopic world we observe.

A (microscopic) particle suspended in a fluid (like dust in air or pollen in water) undergoes apparently random motion due to the many, constant, random collisions of the suspended particle with the fast-moving molecules of the surrounding medium.

We don’t see larger particles jiggle because they are too massive and an collisions aren’t strong enough to noticeable move them.

Although the number of collisions is enormous — on the order of \(10^{16}\) to \(10^{20}\) per second — the forces from these collisions are not perfectly balanced in time or space at every instant. This imbalance creates a net movement, resulting in the irregular, zigzag trajectory known as Brownian motion.

Mathematically, Brownian motion is a continuous stochastic process whose changes over any time interval are drawn from a normal distribution. It is the limit of a random walk as step size and time step go to zero. (The normal distribution is especially useful here because of its favourable analytical properties: it’s stable under addition, has finite variance, and leads to elegant mathematical descriptions of random processes.)

Random walks#

A random walk is a stochastic process in which a particle or agent moves step by step in random directions, with each step independent of the previous ones.

In a random walk as long as the steps are independent, identically distributed (i.i.d.), with finite variance and zero mean then a particle’s expected displacement remains zero, but its typical distance from the origin grows like \(\sqrt{n}\).

In a symmetric random walk:

Each step is equally likely to go left or right (or in random directions in 2D/3D).
So for every path that ends at \(x = +a\), there’s another equally likely path that ends at \(x = -a\).
When you average all those displacements, the positive and negative values cancel out:

\[ \langle x \rangle = 0 \]

But the particle is almost never actually at the origin. It keeps moving — just in random directions.

What grows over time is the spread of the possible positions (the variance), or more precisely:

\[ \langle x^2 \rangle = \text{mean square displacement} \sim n \]

Taking the square root gives the typical distance:

\[ \sqrt{\langle x^2 \rangle} \sim \sqrt{n} \]

So:

The average position stays centred at zero
But the cloud of possible positions grows wider — and the particle wanders further from the origin as time goes on.

Analogy: Imagine a drunk person starting at a lamp post. They take one random step every second. After a minute, they’re just as likely to be 5 metres east as 5 metres west — so on average, they’re still “at the lamp post” in expectation. But in reality, they’re probably somewhere around 7–8 metres away, not still standing at the base.

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.collections import LineCollection
from matplotlib.colors import Normalize
plt.rcParams['text.usetex'] = True

def simulate_fixed_step_brownian(num_steps, step_length):
    angles = np.random.uniform(0, 2*np.pi, size=num_steps)
    displacements = np.column_stack((np.cos(angles), np.sin(angles))) * step_length
    positions = np.cumsum(displacements, axis=0)
    return positions

# Parameters
num_steps = 10000
step_length = 1.0  #alternative is for non-fixed step size e.g. np.random.normal(0, sqrt(2D*dt))
positions = simulate_fixed_step_brownian(num_steps, step_length)

# Create segments between consecutive points
def get_segments(data):
    return np.array([data[i:i+2] for i in range(len(data)-1)])

# Animation setup
fig, ax = plt.subplots()
ax.set_title("Brownian Motion (Fixed Step Size)")
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")

# Centre view on median of path
centre = np.median(positions, axis=0)
half_window = 50  # Adjust based on step length and number of steps
ax.set_xlim(centre[0] - half_window, centre[0] + half_window)
ax.set_ylim(centre[1] - half_window, centre[1] + half_window)

# Prepare the LineCollection
segments = get_segments(positions)
norm = Normalize(0, len(segments))
cmap = plt.cm.Greys
colors = cmap(norm(np.arange(len(segments))))
line_collection = LineCollection([], linewidth=2)
ax.add_collection(line_collection)

plt.close()

# Frame selection
num_frames = 100
frames = np.linspace(2, num_steps - 1, num_frames, dtype=int)

def init():
    line_collection.set_segments([])
    return line_collection,

def update(frame):
    current_segments = get_segments(positions[:frame])
    current_colors = cmap(norm(np.arange(len(current_segments))))
    line_collection.set_segments(current_segments)
    line_collection.set_color(current_colors)
    return line_collection,

animation = FuncAnimation(fig, update, frames=frames, init_func=init, blit=True)

from IPython.display import HTML
HTML(animation.to_jshtml())

With many particles performing random walks like this the particles will slowly spread outwards at a certain rate. This is diffusion.

Eventually they will fill the space they have to occupy and be evenly mixed, i.e. diffusion should lead to uniformity.

Intuitively we know this happens. e.g. you smell dinner before you see it and wearing masks in the midst of a pandemic is good.

The equation Einstein derived to describe Brownian motion is:

\[ \langle x^2 \rangle = 2Dt, \]

where:

\(\langle x^2 \rangle\) is the mean square displacement of the particle (in one spatial dimension)
\(D\) is the diffusion coefficient
\(t\) is the time elapsed
The 2 comes from how variance accumulates in time due to random steps in both directions.

Note that the mean square displacement formula for diffusion tells you how far, on average, a diffusing particle ends up from where it started. It is true for the ensemble average, not a single trajectory.

The diffusion coefficient \(D\) captures the cumulative effects of random molecular collisions at the microscopic level, things like:

The average speed and mass of molecules
The temperature of the fluid, \(T\)
The viscosity of the medium, \(\eta\)
The radius of the diffusing particle, \(r\)

Hence, it can be related to other physical quantities through the Stokes-Einstein equation:

\[ D = \frac{k_B T}{6 \pi \eta r}, \]

where:

\(k_B\) is the Boltzmann constant

By combining this with his earlier result \(\langle x^2 \rangle = 2Dt\), Einstein provided a way to calculate \(D\) (and thus probe the microscopic properties of matter) using macroscopic observations. i.e. the microscopic randomness that we can’t see (molecular collisions wrapped up in \(D\)) result in a macroscopic predictability that can be observed and measured (average behaviour of particle spreading out with time).

This was groundbreaking: in 1905, the existence of atoms and molecules was still debated. Einstein’s work on Brownian motion offered one of the first direct experimental confirmations of atomic theory. By measuring how visible particles diffused through a fluid, he could estimate the size of atoms and even calculate Avogadro’s number. This work helped complete the statistical mechanics programme begun by Maxwell, Boltzmann, and Gibbs — demonstrating how macroscopic phenomena (like diffusion) arise from the statistical behaviour of microscopic particles.

../../_images/c37aa207944ddffbd2f9d6a391e952c2f3dd0fccf985e67b664493d1f0c69ac1.png

Returning to the example we are considering of the two morphogens…

Let the diffusion of \(B\) (a single particle depicted with white-blue-purple), \(D_B\) to be twice that of \(A\) (a single particle depicted with yellow-orange-brown), \(D_A\).

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

# Parameters
n = 2000
diffusion_rate_1 = 0.1  # Particle A (slower)
diffusion_rate_2 = 0.2  # Particle B (faster)

# Simulate Brownian motion
x_1 = np.cumsum(np.sqrt(2 * diffusion_rate_1) * np.random.randn(n))
y_1 = np.cumsum(np.sqrt(2 * diffusion_rate_1) * np.random.randn(n))
x_2 = np.cumsum(np.sqrt(2 * diffusion_rate_2) * np.random.randn(n))
y_2 = np.cumsum(np.sqrt(2 * diffusion_rate_2) * np.random.randn(n))

# Smooth trails (interpolation)
k = 10
x_1B = np.interp(np.arange(n * k), np.arange(n) * k, x_1)
y_1B = np.interp(np.arange(n * k), np.arange(n) * k, y_1)
x_2B = np.interp(np.arange(n * k), np.arange(n) * k, x_2)
y_2B = np.interp(np.arange(n * k), np.arange(n) * k, y_2)

# Compute distance from origin
disp_1 = np.linalg.norm(np.column_stack((x_1, y_1)), axis=1)
disp_2 = np.linalg.norm(np.column_stack((x_2, y_2)), axis=1)

# Reference sqrt curves
sqrt_1 = np.sqrt( 4*diffusion_rate_1 * np.arange(n))
sqrt_2 = np.sqrt( 4*diffusion_rate_2 * np.arange(n))

# Matplotlib colours to match
color_A = '#D55E00'  # reddish-orange (Particle A)
color_B = '#0072B2'  # blue (Particle B)

# Set up side-by-side plots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# --- Left: 2D trajectories ---
sc1 = ax1.scatter(x_1B, y_1B, c=range(n * k), cmap='Oranges', s=1, label="Particle A")
sc2 = ax1.scatter(x_2B, y_2B, c=range(n * k), cmap='Blues', s=1, label="Particle B")

ax1.set_title("2D Trajectories")
ax1.axis('equal')
ax1.set_axis_off()

# Dummy handles for legend
import matplotlib.patches as mpatches
patch1 = mpatches.Patch(color=color_A, label="Particle A (slow)")
patch2 = mpatches.Patch(color=color_B, label="Particle B (fast)")
ax1.legend(handles=[patch1, patch2], loc='upper left', frameon=False)

# --- Right: Distance vs Time ---
ax2.plot(disp_1, label="Particle A", color=color_A)
ax2.plot(sqrt_1, '--', label=r'Ref $\propto \sqrt{t}$ (A)', color=color_A)

ax2.plot(disp_2, label="Particle B", color=color_B)
ax2.plot(sqrt_2, '--', label=r'Ref $\propto \sqrt{t}$ (B)', color=color_B)

ax2.set_xlabel(r"Time step ($t$)")
ax2.set_ylabel(r"Distance from origin ($\|x(t) - x_0\|$)")
# ax2.set_title("Distance from Origin Over Time")
ax2.grid(True)
ax2.legend()

plt.tight_layout()
plt.show()

../../_images/009d2f2edf090d5e30b3e2b2ce58b6a1f21a994f6c3927f0e13ce30d39328b8a.png

Turing instabilities#

Turing examined a chemical system in stable equilibrium, where, in the absence of spatial effects (variations in concentration of particles across space), no net reactions were occurring. He demonstrated that this stability could be disturbed solely by introducing diffusion.

This result is sort of counterintuitive: diffusion is typically associated with smoothing out concentration differences, not amplifying them. But Turing demonstrated that under certain conditions, diffusion can destabilise a homogeneous state and give rise to spatial patterns.

How?… Small random fluctuations in concentration can trigger a reaction–diffusion feedback loop between the interacting activator and inhibitor.

For example, a small, random local spike in activator concentration triggers:

More of the activator to be produced (positive feedback), and
Production of some of the inhibitor.
The activator is slow to diffuse, so it tends to stay concentrated in the region where it was produced. On the other hand, the inhibitor diffuses relatively quickly spreading out faster and suppressing the activator production in surrounding areas, which causes:
The central region to retain high activator levels with the positive feedback loop (making more of itself) continuing locally, as long as it’s not overwhelmed by the returning inhibitor.
Nearby regions are suppressed, creating spatial contrast.

This self-reinforcing interaction between local activation and long-range inhibition causes an initially uniform system to evolve into a stable, patterned structure.

Regulation#

We will also include:

Feed rate, \(f\): controls the supply of \(A\)
Kill rate, \(k\): controls the decay of \(B\) (converted to \(C\) via the chemical reaction: \(B \rightarrow C\), where \(C\) is an inert product of these irreversible reactions)

This is also lumped in with the ‘reaction’ part of the dynamics (because it isn’t diffusion).

Putting it all together#

So what happens in this system?#

How do we find out?

We won’t solve it analytically, in general this is far to complicated. Instead, we simulate it. And that is a kind of solution.

We’ll release particles, allow them to interact according to simple rules, and then observe the resulting behaviour. Patterns — if they emerge — will do so from the dynamics themselves, not from anything we impose.

Assuming we have the right reaction equations, we move on to explore the parameters:

\(k\)
\(f\)
\(D_A\)
\(D_B\)

There is also a length scale for the interaction to occur, possibly some probability for the conversion etc but these are the main parameters.

We need to explore (via parameter sweeps) the many combinations (the parameter space) of these parameters.

4-dimensional parameter space, maybe 3 if we have a relationship between the diffusions. Hopefully we can fix the diffusions using something anchored with physical meaning and then we only need to explore the parameters we have some direct, experimental control over: \(k\) and \(f\).

Not all points in the parameter space produce dynamically interesting patterns.

Turing proposed that the interaction of the morphogens is able to produce six different scenarios upon reaching a stable, constant state:

Case I: Uniform, stationary state: (\(\dot{A}=0, \dot{B}=0\)) unmoving in time and space.
e.g.: Mixing two food colourings in water produces a single, uniform hue.
Case II: Uniform, oscillating state: all cells oscillate together, departing from equilibrium in synchrony.
e.g.: Flowers opening and closing in unison.
Case III: Salt-and-pepper pattern: cells differentiate in a scattered fashion, inhibiting nearby cells from doing the same. e.g.: Spacing of sunflowers in a garden due to local competition for resources.
Case IV: Oscillating salt-and-pepper: the salt-and-pepper pattern itself varies over time.
Case V: Travelling wave: Waves of concentration move through the system (usually requires three or more morphogens).
e.g.: A slinky wave propagating from one end to the other.
Case VI: Stationary wave (Turing pattern): stable, repeating pattern emerges and persists despite ongoing reactions and diffusion.
e.g.: Zebra stripes or fingerprint ridges, where the wavelength depends on reaction and diffusion rates.

Ultimately, the best way to understand how these parameters shape the system’s behaviour is to run a simulation, like this one:

Complexity Explorables

Have you seen the Complexity Explorables yet? Warning: You may lose a big part of your day visiting this website!
Another really nice one
And the ‘original’ here Karl Sims’ had an exhibit in the Museum of Science in Boston that was based on this tool

Extensions#

The reaction-diffusion system can be modulated by all sorts of other factors that will influence the formation of the patterns but require extending (and complicating) the model, e.g.:

Orientation: diffusion can occur faster in one direction than another to give an orientation to the results.
‘Style Map’: the feed and kill rates can vary across the grid to give different patterns in different areas.
Flow: the chemicals can flow across the grid under the influence of an external force to give various dynamic effects.
Size: the scale of the pattern changes when the reaction rate is sped up or slowed down relative to the diffusion rate.
External influences

Returning to why a cheetah has spots#

Many animals display distinct and repeatable patterns—stripes, spots, or swirls—yet they share a remarkably similar set of genes. So how does this visual diversity arise?

The key lies not in each species having a unique genetic “blueprint” for their patterns, but in how biological shape and scale influence the behaviour of how the morphogens spread and interact to stimulate or suppress development.

The beauty of mathematical modelling is that once you're confident your model captures the biology, you can start to play around with it and ask biological questions that are difficult to answer through experiments alone.

— Kit Yates

Their skin is patterned with horizontal stripes made not from chemical concentration directly, but by the spatial arrangement of pigment cells:

Melanophores (black)
Xanthophores (yellow)
Iridophores (silvery-blue)

These cell types interact, migrate, and compete in ways that mimic a reaction-diffusion system. A study by Kit Yates and colleagues used mathematical modelling to understand which biological processes were essential for stripe formation and found that the orientation of the stripes—horizontal rather than vertical—emerged naturally from the geometry of the fish.

Returning to how zebras — or sand dunes — get their stripes#

In modelling, we often use a familiar system to understand another with similar underlying dynamics. For instance, a mechanical mass-spring-damper system can help us understand an electrical RLC circuit because they follow the same mathematical form with both systems governed by second-order linear ODEs and exhibit oscillatory or overdamped behaviour depending on the parameters:

Mechanical system: Mass–Spring–Damper: a mass \(m\) attached to a spring (stiffness \(k\)) and damper (damping coefficient \(c\)) obeys:

\[ m \ddot{x}(t) + c \dot{x}(t) + k x(t) = F(t) \]

\(x(t)\): displacement of the mass
\(\dot{x}(t), \ddot{x}(t)\): velocity and acceleration
\(F(t)\): external force applied

Electrical system: RLC Circuit: an RLC series circuit with resistance \(R\), inductance \(L\), and capacitance \(C\) obeys:

\[ L \ddot{q}(t) + R \dot{q}(t) + \frac{1}{C} q(t) = V(t) \]

\(q(t)\): electric charge
\(\dot{q}(t) = i(t)\): current
\(V(t)\): applied voltage

Analogy mapping:

Mechanical (mass-spring-damper)	Electrical (RLC circuit)
Displacement \(x\)	Charge \(q\)
Velocity \(\dot{x}\)	Current \(i = \dot{q}\)
Force \(F(t)\)	Voltage \(V(t)\)
Mass \(m\)	Inductance \(L\)
Damping coefficient \(c\)	Resistance \(R\)
Spring constant \(k\)	Inverse capacitance \(1/C\)

In the same spirit, if the chemistry behind zebra stripe pigmentation feels abstract, we might instead turn to something more tangible, like sand dunes to provide an alternative lens:

Wind-blown sand forms ripples through self-reinforcing feedback
once a ripple forms, it captures more sand and grows.
At the same time, this growing ripple depletes sand from the wind stream, suppressing the formation of another ripple nearby.
This dynamic interplay of local activation (more sand builds a bigger ripple) and long-range inhibition (downwind suppression) echoes the same principles behind Turing’s chemical model of pattern formation.

Turing-like behaviour isn’t limited to biochemistry. Similar activator-inhibitor dynamics have been proposed in:

Ecology (e.g., predator-prey systems)
Neuroscience (neurons firing or inhibiting neighbours)
Cardiology (electrical wave patterns in the heart)
Sand dunes (ripples from wind as a diffusive force)
Ant colonies (structures formed from arranging corpses or colony trails)
Finger prints (mouse digits were used as a model for understanding the formation of human fingerprints)
Hallucinated visual patterns (grids, spirals, tunnels may arise from neural Turing mechanisms in the visual cortex.)

Turing’s model marked a pivotal shift in biology, introducing a more mechanistic and mathematical approach to development. By treating biological pattern formation as a physical process governed by simple rules and constraints, he laid the groundwork for modelling systems that had previously seemed too messy, nonlinear, or complex to capture mathematically.

Of course, real biology is far more intricate than Turing’s original two-morphogen setup. Developmental processes typically involve dozens to hundreds of interacting compounds, operating within a highly heterogeneous and dynamic environment. Identifying which molecules act as activators or inhibitors remains a significant challenge—many of the relevant kinetic parameters are difficult, if not impossible, to measure directly in living tissue.

Despite the elegance of the theory, no universally accepted molecular Turing mechanism has been identified in natural developmental systems. The patterns we observe—whether in skin pigmentation, limb formation, or cellular arrangement—are likely the result of multiple overlapping processes, including:

Reaction-diffusion dynamics
Mechanical feedback and tissue deformation
Cell growth and migration
Genetic regulation and signalling cascades

and Turing’s model should be seen not as the whole explanation, but as a powerful component within a larger biological system.

I’ve sort of lied to you…#

The reaction-diffusion system we described involves a huge number of particles.

Tracking each individual particle’s movement for many many generations is enormously computationally expensive.

It can be done (with an agent-based modelling framework like what we will see in a few weeks) but it’s ridiculously complicated (which is why statistical physics aggregates all the particles into average/global properties like temperature). More importantly, it isn’t necessary to reproduce the observed behaviour.

The solutions you saw in the Complexity Explorable were actually the result of a ‘coarse-grained’ model where instead of caring about the individual particles we will lump many of them together into cells and track the cell concentration instead.

Remember…#

Fig from: The first paragraph of Turing’s paper on morphogenesis.

The modeller’s purpose is absolutely, definitely not to reproduce the system exactly

Instead, we want:

to capture the essence of a system
for simulations to run “quickly”
for the model to scale well to larger inputs

Is the world discrete or continuous?#

If it is discrete, when is pretending it is continuous acceptable/helpful?

If it is continuous, when is pretending it is discrete acceptable/helpful?

The true nature of reality might encompass both aspects, and our understanding is still evolving as scientific theories develop.

We have theories for both scenarios:

Discrete Mathematics: deals with countable, distinct elements e.g. integers, graphs, and combinatorial structures
- Quantum Physics: discreteness at fundamental level e.g. energy levels are quantised, particles exhibit both wave-like (continuous) and particle-like (discrete) properties. Space and time may have a discrete structure at the Planck scale?
Continuous Mathematics: deals with uncountable sets and smooth changes e.g. real numbers, calculus, and DEs
- Classical Mechanics and Electromagnetism: treat space, time, and fields as continuous
- General Relativity: describes gravity in terms of the continuous curvature of spacetime

The modelling tools and frameworks you choose depend on several factors:

Scale: continuous models often work well at macroscopic scales, while discrete models are better suited to microscopic or agent-level dynamics.
Nature of the phenomenon: some systems are inherently discrete (e.g. individuals, events), others continuous (e.g. temperature, concentration).
Level of detail required: finer-grained models capture more nuance but may obscure general patterns or become intractable.
Computational constraints and resources available: in practice, we are limited by resolution and memory (and in the end, our predictions are only as accurate as the pixels we can compute).

Is Brownian motion discrete or continuous?#

Aspect	Nature	Explanation
Particle	Discrete	You’re tracking the motion of an individual microscopic particle (e.g. a pollen grain), which is a discrete object.
Space	Continuous	The particle’s position \(x(t)\) is represented as a real-valued variable in continuous space.
Time	Continuous (in theory)	Although early explanations used time steps, the modern model treats time as continuous.
Motion	Continuous but non-smooth	The motion is continuous in time but highly irregular — it’s modelled as a Wiener process, which is continuous everywhere but differentiable nowhere.

Modelling Brownian motion in a discrete world#

Brownian motion — the jittery movement of a microscopic particle suspended in fluid — can be understood as the result of many small, random steps.

This is by a random walk.

The random walk is a Markov process with each step depending only on the current state, not the past (memoryless).

Imagine a particle that every \(\Delta t\), jumps left or right by a distance \(\Delta x\), chosen at random… this offers a discrete, intuitive model of particle motion.

While it’s theoretically possible to model every particle’s path using Newton’s laws combined with random forces, doing so is often unnecessary and computationally expensive. In many cases (especially when we’re interested in average behaviour or large-scale patterns) it’s far more efficient to use continuum-scale models. Turing, for example, wasn’t concerned with individual particle trajectories but with evolving patterns of concentration.

Modelling diffusion in a continuous world#

The alternative is to aggregate particles.

If we do this, then the spatio-temporal dynamics can be modelled with partial differential equations in a continuous field.

Diffusion arises from the random motion and collisions of particles at the microscopic level. This causes particles to spread from regions of high concentration to regions of lower concentration, gradually smoothing out the concentration field. We can describe this behaviour using a flux that is proportional to the concentration gradient, and track how the concentration changes over time as a result of this particle movement.

Let \(P(x, t)\) be the probability the particle is at position \(x\) at time \(t\). Over time, this distribution spreads out. Thanks to the Central Limit Theorem, the sum of many such random steps converges to a normal distribution. This explains why displacements in Brownian motion are Gaussian-distributed, even if step sizes are fixed.

Einstein formalised this in 1905, showing that the mean square displacement grows linearly with time:

\[ \langle x^2 \rangle = 2Dt, \]

where:

\(\langle x^2 \rangle\) is the mean square displacement
\(D\) is the diffusion coefficient

The key relationship linking space and time in a discrete random walk is:

\[ D = \frac{(\Delta x)^2}{2\Delta t} \]

This defines the diffusion coefficient in terms of the size and frequency of the steps.

Even though particles are discrete and their motion is random, smooth and predictable laws emerge when we look at the behaviour over time or across many particles.

If we let the step size and time step become infinitesimally small — while keeping \(D\) fixed — we get a smooth, continuous limit:

\[ \frac{\partial P(x,t)}{\partial t} = D \frac{\partial^2 P(x,t)}{\partial x^2} \]

This is the diffusion equation, which takes the same mathematical form whether we are modelling:

\(P(x,t)\): the probability density of a single particle’s position, or
\(C(x,t)\): the concentration of many particles over space and time.

The power of statistical mechanics is that deterministic equations like the diffusion equation can describe the statistical behaviour of stochastic systems. Einstein’s insight was in linking how those deterministic models of the macroscopic behaviour arise from the randomness of the underlying system.

e.g., we can consider a flux of particles that is proportional to the concentration gradient and describe how that concentration changes with time as particles diffuse from regions of high concentration to regions of low concentration.

Gif from Wiki: Molecular diffusion from a microscopic and macroscopic point of view

Fick’s Laws#

So far, we’ve approached diffusion from a microscopic, probabilistic point of view — showing how the random motion of particles gives rise to a smooth, predictable evolution over time. But there’s another way to arrive at the diffusion equation — one that starts from physical intuition and macroscopic quantities like flux and concentration gradients.

This perspective, due to Fick, is grounded in mass conservation and is widely used in fields like biology, chemistry, and engineering.

Note: If you’re already comfortable with the diffusion equation from the Einstein/random walk view, feel free to skim this section. But Fick’s laws offer valuable intuition and will deepen your understanding of how diffusion is modelled and applied across disciplines.

Fick’s First Law (steady state diffusion): relates the diffusive flux to the concentration gradient. It states that the flux of particles, (J), is proportional to the negative gradient of the concentration, (C): \( J = -D \frac{dC}{dx}\)

\(J\) is the diffusion flux (amount of substance per unit area per unit time),
\(D\) is the diffusion coefficient (a measure of how quickly particles diffuse),
\(\frac{dC}{dx}\) is the concentration gradient,
minus sign says that diffusion goes from high to low concentration

Fick’s Second Law (time-dependent diffusion): describes how the concentration of particles changes with time as proportional to the local curvature of the concentration gradient. It is derived from the first law and the principle of mass conservation: \(\frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2}\)

\(\frac{\partial C}{\partial t}\) is the rate of change of concentration with time,
\(\frac{\partial^2 C}{\partial x^2}\) is the second spatial derivative of the concentration.

The diffusion equation (aka heat equation) can then be derived from Fick’s laws combining how flux responds to gradients (1st law) and how that flux changes concentration over time (2nd law):

\[ \frac{\partial C}{\partial t} = -\frac{\partial}{\partial x} \left( -D \frac{\partial C}{\partial x} \right) = D \frac{\partial^2 C}{\partial x^2} \]

In two or more dimensions we generalise the derivatives with the gradient and Laplacian operators: \(J=−D\nabla C\)

\(\frac {\partial C }{\partial t}=D\nabla^2 C. \)

Recall that:

The gradient operator (a vector) is defined as:

\[ \nabla =\Big( \frac{\partial }{\partial x_1}, \frac{\partial }{\partial x_2}, \ldots, \frac{\partial }{\partial x_n}\Big). \]

The Laplacian, \(\nabla^2\) (or \(\Delta\)), operator is defined with the second partial derivatives of a function with respect to its spatial variables:

\[ \nabla^2 f=\frac{\partial^2 f}{\partial^2 x_1} + \frac{\partial^2 f}{\partial^2 x_2} + \ldots + \frac{\partial^2 f}{\partial^2 x_n}. \]

This is a generalisation of a second derivative of a function (second-order because it is the divergence of a gradient field) and can be interpreted in a similar way, i.e.

how the function changes ‘on average’ as you move away from a given point (positive/negative Laplacian values denote an average increase/decrease respectively)
it captures the spatial variation and curvature of a function by quantifying the rate of change of a scalar field

Check your foundational understanding:
What exactly is the gradient? What exactly is the Laplacian?

The gradient (\(\nabla f\)) is a vector that points in the direction of steepest increase of a scalar field and gives the rate of change in that direction. The Laplacian (\(\nabla^2 f\)) is a scalar that measures the net rate of flow in or out of a point — i.e. the sum of the second derivatives. It tells you if the point is a peak, valley, or flat relative to its surroundings.

In short:

Gradient: “Which way is up?”
Laplacian: “Am I higher or lower than my neighbours?”

Transport equation#

The transport equation is a more general form that describes the movement of substances (such as particles, heat, or momentum) in a medium. It incorporates additional factors such as convection, advection, sources/sinks and reaction, etc and therefore is the umbrella under which Fick’s laws, the diffusion equation, and Turing’s reaction–diffusion model all live.

\[ \frac{\partial u}{\partial t} + \nabla \cdot (\mathbf{v} u) = D \nabla^2 u + R(u) \]

\(u(\mathbf{x}, t)\): the scalar quantity being transported (e.g. concentration)
\(\mathbf{v}\): velocity field
\(D\): diffusion coefficient
\(R(u)\): source, sink or interaction term
\(\nabla \cdot (\mathbf{v} u)\): advection (transport by flow)
\(D \nabla^2 u\): diffusion (spreading out due to random motion)

Special cases of the transport equation:

\(\frac{\partial u}{\partial t} + \mathbf{v} \cdot \nabla u = 0\): advection only, transport by flow, no diffusion or reaction
\(\displaystyle \frac{\partial u}{\partial t} = D \nabla^2 u\): diffusion only, passive spreading due to diffusion (Fick’s Law)
\(\frac{\partial u}{\partial t} = D \nabla^2 u + R(u)\): reaction–diffusion system, no flow

Return to the reaction-diffusion system#

It is the interplay between the reaction and diffusion that creates emergent patterns.

Reaction-diffusion equations are a class of models that separate the non-spatial and spatial dynamics by describing them via reaction terms and diffusion terms respectively:

\[ \frac{\partial U}{\partial t}=R(U) + D \nabla^2 U. \]

reaction terms, \(R(U)\), describe the local dynamics and how different chemical concentrations are interacting, e.g. reaction rates, production or consumption rates, sources, sinks, decay etc.
diffusion terms, \(D \nabla^2 U\), are strictly limited to the Laplacian of the state variable.

A reaction-diffusion system then consists of multiple coupled reaction-diffusion equations:

\[ \frac{\partial U}{\partial t}=R_U(U,V) + D_U \nabla^2 U \]

\[ \frac{\partial V}{\partial t}=R_V(U,V) + D_V \nabla^2 V \]

with more concentrations and further coupling defined similarly if needed.

Recall that in the Gray Scott system, the activator, \(A\), and inhibitor, \(B\) react as:

\[ A + 2B\rightarrow 3B \]

\(U\) and \(V\) are the concentrations of \(A\) and \(B\) respectively.

“Solving” the PDEs#

Solving a system of reaction–diffusion equations analytically is often impossible due to their nonlinear and spatially extended nature. Instead, we simulate these systems numerically to explore how patterns emerge over time from the interplay of local reactions and diffusion across space.

Our goal is to construct a “coarse-grained” model that enables us to observe the emergence of Turing patterns without the computational overhead of monitoring individual particles.

Effective simulation of complex systems requires moving fluidly between discrete and continuous models:

The chemical reactions and diffusion-driven pattern formation is happening at the particle level (discrete).

But… Spatio-temporal dynamics is typically described via particle concentrations with PDEs (continuous).

But… We need to actually simulate this (discrete again!).

Diffusion in a coarse-grained (discretised) world…#

In a discrete setting, we represent space and time as divided into localised, distinct units:

Discrete space: the domain is divided into a grid, and particles are aggregated within cells to represent local concentrations.
Discrete time: the system is evolved in fixed time steps, allowing us to iteratively update concentrations.

To model diffusion in this coarse-grained world, we approximate continuous mathematical operators — such as the gradient and divergence — using finite differences. These discrete operators describe how quantities like concentration change locally across space and time.

This doesn’t require full calculus: partial differential equations can be translated directly into code using simple update rules. e.g., the concentration fields \(U\) and \(V\) evolve over time as:

\[ U \rightarrow U + D_U \, \Delta U \]

\[ V \rightarrow V + D_V \, \Delta V \]

Here, \(\Delta U\) and \(\Delta V\) represent the change in concentration due to diffusion, often approximated numerically using the discrete Laplacian, denoted \(\nabla^2 U\) and \(\nabla^2 V\).

Before introducing formal operators like the Laplacian, we can explore diffusion using an intuitive, rule-based approach.

Consider a simple scenario:

We represent space as a frid, where each cell holds a concentration of particles of type \(A\), denoted by \(U\).
Initially, only the centre cell contains particles, with \(U = 1\). All other cells have \(U = 0\).
At each time step, we assume that 20% of the particles in a cell diffuse equally to their four von Neumann neighbours (North, South, East, West).
Total concentration is conserved across the grid.

This rule gives us a way to simulate diffusion without using derivatives. By updating the grid step by step, we can observe how concentration spreads out over time — capturing the essential behaviour of diffusion with simple arithmetic.

Things (particles/concentration/\(U\)) are moving quickly here. We can slow it down by multiplying by a factor (the diffusion rate, \(D_U\)). For example, setting \(D_U=0.1\) will mean that only 10% of the central cell’s concentration will be committed to contribute to the spread at a given time.

We can visualise the numbers (representing concentration) in these cells as colour. This is an important abstraction as the domain gets large (or the cell gridding becomes fine)

As time evolves, we can animate the grid to see how the concentrations are changing over time.

This depiction focusses on a single particle concentration. If there are more than one type of particle then these can be coarse grained separately such that a single cell contains information about each concentration.

You can then visualise the two concentrations separately or they could be combined in some sensible way e.g. as a proportion \(\frac{U}{U+V}\).

Practice:
Simulate diffusion on a 2D grid by applying the following rule: at each time step, assume that 25% of the concentration at each grid cell spreads equally to its four von Neumann neighbours (up, down, left, right).

Evolve the system over several time steps and observe how the concentration field smooths out. Do this by hand and then, implement the same update in Python and verify that your simulation matches your manual result.

../../_images/eb0cbd124271892a46e7d9da1a9f575b0c7e0df25f7d9928390893a967983815.png

A note on choosing colour schemes#

I did my PhD (and continue to work) with someone who is red-green colour blind. It means he can’t see the number sitting inside this circle:

Whilst I am marking your projects and can see a full rainbow, please ensure you are as inclusive as possible with your writing and your figures.

See here and here for a brief guide on selecting colour schemes

Numerical approximation of the gradient using the finite-difference method#

Our simple rule-based update gives us a good intuitive feel for how diffusion spreads concentration across space.

To build more flexible, accurate, and scalable models, we use discrete approximations of key operators — the gradient and Laplacian — to describe how quantities change across space and time on a grid.

To approximate the gradient of a scalar field \(U\) on a grid, we use the finite-difference method, which builds on the definition of the derivative:

\[ \frac{df}{dx} = \lim_{h \to 0} \frac{f(x + h) - f(x)}{h} \]

In a discrete setting with uniform grid spacing \(h = 1\), we approximate the partial derivatives using central differences, which are more accurate and spatially symmetric.

At grid cell \((i, j)\), the partial derivatives are:

\[ \frac{\partial U}{\partial x}(i,j) \approx \frac{U(i+1,j) - U(i-1,j)}{2} \]

\[ \frac{\partial U}{\partial y}(i,j) \approx \frac{U(i,j+1) - U(i,j-1)}{2} \]

Thus, the gradient vector at \((i,j)\) is approximated by:

\[ \nabla U(i,j) \approx \left( \frac{U(i+1,j) - U(i-1,j)}{2},\; \frac{U(i,j+1) - U(i,j-1)}{2} \right) \]

This gives the rate of change of \(U\) in the \(x\)- and \(y\)-directions respectively, based on values from both sides of the point.

The central difference scheme ensures no directional bias and is second-order accurate. Forward or backward differences can be used instead when information is only available on one side, e.g. near boundaries or in time-dependent problems where causality matters. These one-sided approximations are simpler but typically less accurate (first-order).

Numerical approximation of the Laplacian with the finite-difference method#

To approximate the Laplacian \(\nabla^2 U\) of a scalar field \(U\) on a grid, we use the finite-difference method to compute second-order partial derivatives in both spatial directions.

For numerical symmetry, we approximate the second derivative using values on both sides of the point of interest. This ensures the update depends equally on neighbours in the positive and negative directions — which is important in space, where no direction is preferred. Note: Unlike space, time has a preferred direction (the “arrow of time”), so symmetry isn’t required when approximating time derivatives.

In one dimension, the second derivative is approximated by:

\[ \frac{\partial^2 U}{\partial x^2}(i) \approx U(i+1) - 2U(i) + U(i-1) \]

In two dimensions, we compute the second derivatives in both directions:

\[ \frac{\partial^2 U}{\partial x^2}(i,j) \approx U(i+1,j) + U(i-1,j) - 2U(i,j) \]

\[ \frac{\partial^2 U}{\partial y^2}(i,j) \approx U(i,j+1) + U(i,j-1) - 2U(i,j) \]

Adding these gives the discrete Laplacian:

\[ \nabla^2 U(i,j) \approx U(i+1,j) + U(i-1,j) + U(i,j+1) + U(i,j-1) - 4U(i,j) \]

This expression computes the difference between the value at the centre cell and the average of its North, South, East, and West neighbours — known as the von Neumann neighbourhood or the five-point stencil. The term \(−4U(i,j)\) ensures the centre is weighted to balance the four surrounding values symmetrically. i.e. it gives a measure of how different the value at \((i,j)\) is from the average of its four nearest neighbours — a key driver of diffusion-like behaviour.

This operation can be expressed as a convolution with the following kernel (or stencil matrix):

\[\begin{split} L=\left[ \begin{array}{ccc} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{array} \right] \end{split}\]

This discrete Laplacian matrix serves the same purpose as the continuous Laplacian operator — quantifying local spatial variation — and is commonly used in modelling diffusion.

Theoretical Question:
Why is this diffusion?

Because it’s the difference between the average of nearby grid cells and the central cell of interest (a local averaging at every point in the domain). This causes the cell to become more like their neighbours, hence simulating diffusion

Practice:
Now implement diffusion using a discrete approximation of the Laplacian operator:

\[ \nabla^2 u(i, j) \approx u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1) - 4u(i,j) \]

Use this to evolve the system over time with diffusion coefficient \(D_U=0.25\)

First, compute one time step by hand for a small grid.
Then, implement the update in Python using a manual nested loop approach.
Next, repeat the same update using a convolution with the five-point stencil.
Compare the results of both implementations. Do they produce the same outcome?
How do these compare to the earlier rule-based diffusion method? What advantages do you notice in using the Laplacian formulation, and what are the trade-offs between manual and convolution-based approaches?

Regulation#

In addition to diffusion, chemical reactions regulate the local concentrations of particles by adding, removing, or transforming them. These processes happen can be incorporated directly into the model as additional update rules at each time step.

feed:#

\(A\) particles are fed into the system. Let’s assume that particles are added to each cell, at each time step.

Adding particles at a constant rate \(f\) to the concentration may cause \(U\) to exceed our maximum value.

If we instead add particles at a rate \(f(1-U)\) then the amount added depends on the current concentration of the particular grid cell and is limited to the maximum normalised value of 1.

\[ U(i,j)\rightarrow U(i,j) + f(1-U(i,j)) \]

We also remove particles \(B\) (and \(C\)) at a rate proportional to their current concentration of the particular grid cell to keep the total density of particles fixed:

\[ V(i,j)\rightarrow V(i,j) - fV(i,j) \]

kill:#

\(B\) particles are converted to inert \(C\) particles at a rate, \(k\), proportional to their current concentration at the particular grid cell.

\[ V(i,j) \rightarrow V(i,j) - kV(i,j) \]

Reaction#

The likelihood of a particular reaction occurring depends on how frequently particles collide within a given grid cell — which in turn depends on their local concentrations.

Consider:

\[ A + 2B\rightarrow 3B \]

Since two \(B\) particles are required, the reaction is more sensitive to the concentration of \(B\) than \(A\). For example, when \(V\) is small, the probability of finding two \(B\) particles in the same cell drops off quickly. To reflect this, we model the reaction rate as being proportional to:

\[ UV^2 \]

That is, the rate increases linearly with \(U\) and quadratically with \(V\). Once such a collision occurs, the reaction may only proceed at some probability \(r \in [0,1]\). For simplicity, we set \(r = 1\).

The resulting updates to the concentrations are:

\[ U \rightarrow U - UV^2 \]

\[ V \rightarrow V + UV^2 \]

(I’ve omitted the \((i,j)\) subscripts here for clarity, but note that these updates are applied locally at each grid cell.)

Putting it all together:#

The reaction + diffusion of particles was modelled with PDEs:

\[ \frac{\partial U}{\partial t}= D_U \nabla^2 U + f(1-U)-UV^2 \]

\[ \frac{\partial V}{\partial t}=D_V \nabla^2 V - (k+f)V +UV^2 \]

that were discretised for numerical simulation:

\[ U \rightarrow U + D_U\Delta U + f(1-U) - UV^2 \]

\[ V \rightarrow V + D_V\Delta V - (k+f)V+ UV^2 \]

This is the Gray-Scott model - one of the most well studied reaction-diffusion systems.

Don’t let their simplicity fool you - the interplay between reaction and diffusion of ‘morphogens’ can give rise to spatio-temporal dynamics that is highly complex and diverse.

Simulating#

Models like this can’t be solved analytically.

But that’s ok because we don’t really want to solve it.

i.e. we don’t really care where the stripes on a zebra are, only that for a particular set of parameters and conditions stripes occur.

There are some beautiful simulators out there to play with. See for example:

But we would still like to create our own… You’ll do this week’s Workshop.

So is the world discrete or continuous?#

We began by asking whether the world is fundamentally continuous or discrete — and through the lens of Brownian motion, we’ve explored how both perspectives can be valid, depending on how we model and at what scale we observe.

Modelling Brownian motion with a random walk#

This is a microscopic and stochastic description of the motion of an individual particle as a series of discrete, random steps in space and time.

Aspect	Nature	Explanation
Time	Discrete	Einstein modelled Brownian motion as a sequence of random jumps in small time intervals.
Steps	Discrete	In a random walk, the particle takes discrete steps in space at each discrete time interval.
Transition	Probabilistic	Each step is random — often modelled as drawn from a normal distribution with zero mean, leading to no net drift over time.
Limit	Continuous stochastic process	As the time step and the step length scales appropriately, the random walk converges to Brownian motion, which is the basis of the diffusion equation.

Modelling Brownian motion with the diffusion equation#

This is a macroscopic and deterministic description of how a distribution of particles evolves smoothly over time in continuous space.

Aspect	Nature	Explanation
Time	Continuous	Time is treated as a continuous variable, allowing smooth evolution of the system over any time scale.
Space	Continuous	Space is modelled as a continuous domain; the concentration of particles is defined at every point in space.
State	Density field	Instead of tracking individual particles, the model describes the concentration of many particles as a scalar field \(u(x,t)\).
Transition	Deterministic	The diffusion equation (a partial differential equation) describes how the concentration evolves over time based on physical laws, such as Fick’s law.
Limit	Emerges from stochastic behaviour	The diffusion equation arises as the continuum limit of a large number of independent random walks, capturing the average behaviour of many particles undergoing Brownian motion.

Week 3: Reaction-diffusion systems

Contents

Week 3: Reaction-diffusion systems#

How the cheetah got his spots#

Patterns#

What exactly is a pattern?#

How do patterns form?#

Alan Mathison Turing#

The chemical basis of morphogenesis#

Reaction-diffusion systems#

Reaction#

Diffusion#

Random walks#

Turing instabilities#

Regulation#

Putting it all together#

So what happens in this system?#

Extensions#

Returning to why a cheetah has spots#

Returning to how zebras — or sand dunes — get their stripes#

I’ve sort of lied to you…#

Remember…#

Is the world discrete or continuous?#

Is Brownian motion discrete or continuous?#

Modelling Brownian motion in a discrete world#

Modelling diffusion in a continuous world#

Fick’s Laws#

Transport equation#

Return to the reaction-diffusion system#

“Solving” the PDEs#

Diffusion in a coarse-grained (discretised) world…#

A note on choosing colour schemes#

Numerical approximation of the gradient using the finite-difference method#

Numerical approximation of the Laplacian with the finite-difference method#

Regulation#

feed:#

kill:#

Reaction#

Putting it all together:#

Simulating#

So is the world discrete or continuous?#

Modelling Brownian motion with a random walk#

Modelling Brownian motion with the diffusion equation#